home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / toolbox / house.r < prev    next >
Text File  |  1994-06-26  |  2KB  |  61 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Create a Householder matrix.
  4.  
  5. // Syntax:    house ( X )
  6.  
  7. //    If HOUSE(x), which returns a list containing elements `v' and
  8. //    `beta',  then H = EYE - beta*v*v' is a Householder matrix such
  9. //    that Hx = -sign(x(1))*norm(x)*e_1. 
  10. //    NB: If x = 0 then v = 0, beta = 1 is returned.
  11. //            x can be real or complex.
  12. //            sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).
  13.  
  14. //    Theory: (textbook references Golub & Van Loan 1989, 38-43;
  15. //             Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
  16. //      Hx = y: (I - beta*v*v')x = -s*e_1.
  17. //      Must have |s| = norm(x), v = x+s*e_1, and
  18. //      x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
  19. //      So take s = sign(x(1))*norm(x) (which avoids cancellation).
  20. //      v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
  21. //          = 2*norm(x)*(norm(x) + |x(1)|).
  22.  
  23. //    References:
  24. //        G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  25. //           Johns Hopkins University Press, Baltimore, Maryland, 1989.
  26. //        G.W. Stewart, Introduction to Matrix Computations, Academic Press,
  27. //           New York, 1973,
  28. //        J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
  29. //           Press, 1965.
  30.  
  31. //    This file is a translation of house.m from version 2.0 of
  32. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  33. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  34.  
  35. //-------------------------------------------------------------------//
  36.  
  37. house = function ( x )
  38. {
  39.   local (m, n, s, v, beta)
  40.  
  41.   m = x.nr; n = x.nc;
  42.   if (n > 1) { error ("Argument must be a column vector."); }
  43.  
  44.   s = norm(x,"2") * (sign(x[1]) + (x[1]==0)); // Modification for sign(0)=1.
  45.   v = x;
  46.   if (s == 0)    // Quit if x is the zero vector.
  47.   {
  48.     beta = 1; 
  49.     return << beta = beta ; v = v >>;
  50.   }
  51.  
  52.   v[1] = v[1] + s;
  53.   beta = 1/(s'*v[1]);        // NB the conjugated s.
  54.  
  55.   // beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
  56.   // But beta as above can be non-real (due to rounding) only 
  57.   // when x is complex.
  58.  
  59.   return << beta = beta ; v = v >>;
  60. };
  61.